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Oh; 
"^ , We study, through the diffusion Monte Carlo method, a spin one-half fermion fluid, in 



the three dimensional Euclidean space, at zero temperature. The point particles, im- 
mersed in a uniform "neutralizing" background, interact with a pair-potential which 
can be continuously changed from zero to the Coulomb potential depending on a pa- 



u 

c/3 ■ rameter jj,. We determine the radial distribution functions of the system for various 



c^ . values of density, /i, and polarization. We discuss about the importance, in a com- 

I I puter experiment, of the choice of suitable estimators to measure a physical quantity. 

Ch ' The radial distribution function is determined through the usual histrogram estima- 
O ■ 

, ^, ! tor and through an estimator determined via the use of the Hellmann and Feynman 

theorem. In a diffusion Monte Carlo simulation the latter route introduces a new bias 

00 ■ to the measure of the radial distribution function due to the choice of the auxiliary 

a^ 

CN ' function. This bias is independent from the usual one due to the choice of the trial 
l> 



wave function. A brief account of the results from this study were presented in a 
recent communication [R. Fantoni, Solid state Communications, 159, 106 (2013)]. 
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I. INTRODUCTION 

The Jellium model is a system of pointwise electrons of charge e and number density 
n in the three dimensional Euclidean space filled with an uniform neutralizing background 
of charge density —en. The zero temperature, ground-sate, properties of the statistical 
mechanical system thus depends just on the electronic density n or the Wigner-Seitz radius 
Tg = (3/47m)^/^/ao where Uq is Bohr radius. The model can be used for example as a first 
approximation to describe free electrons in metallic elementsi (2 < r^ < 4) or a white dwar^ 
(r, ~ 0.01). 

When an impurity of charge q is added to the system, the screening cloud of electrons 
will experience the Friedel oscillations. In the Thomas-Fermi description of the static 
screening an electric potential qvH{r) (the Hartree potential) is created by the impurity 
and by the redistribution of the electronic charge n(r) — n. It obeys the Poisson equation 
geV^f_H'(r) = 47re[—q6{r) — en(r) +en] and the equilibrium condition on the electrochemical 
potential, /ic(n(r)) + qevnir) = constant. An analytic solution can be obtained for |g| ^ 1, 
when we find n{r) — n ^ —qevHir)dn/dfic by expansion of fic around the homogeneous 



state. Assuming dn/d^c is positive and with the definition ks = ^/ATre^dn/d^c the Poisson 
equation yields 



e ^"'^ 



VH{r) = . (1) 

r 

It is clear from this result that the quantity l/kg measures the distance over which the self 
consistent potential associated with the impurity penetrates into the electron gas. Thus, 1/ks 
has the meaning of a screening length. The Thomas-Fermi value of the screening length is 
obtained by replacing the thermodynamic quantity dn/d^c by its value for non-interacting 
fermions, using for /i^ the Fermi energy. Clearly we have that f/f(r) -^ 1/r as l/kg — t- cxd 
and vh{j) — ;• as 1/ks — > 0. Also vh is short ranged. 

It is important to study the ground-state properties of a model of point fermions of spin 
one-half interacting with a bare pair-potential v^{r) which can be continuously changed from 
zero (/i — !■ 0, ideal gas) to the Coulomb potential (/i — )■ oo, Jellium model) depending on a 
parameter /i. And we chose the following functional form 

, , erf(ur) , , 

^m(^) = -^ ■ (2) 



Still the fluid is immersed in a static uniform background of continuously distributed point 
particles which interact with the particles of the fluid with the same pair-potential but of 
opposite sign. 

A major challenge in the Kohn-Sham scheme of Density Functional Theory is to de- 
vise approximations to the exchange-correlation functional that accurately describes near- 
degeneracy or long-range correlation effects such as van der Waals forces. Among recent 
progresses to circumvent this problem, we mention "range-separated" density functional 
schemes which combine the Kohn-Sham formalism with either random-phase approximation^ 
or multideterminantal approaches^. Such schemes require a local density functional for par- 
ticles interacting via modified potentials defined in terms of a suitable parameter /i which 
either softens the core or suppresses the long-range tail. Further insight into electronic cor- 
relations in molecules and materials can be gained through the analysis of the on-top pair 
correlation function^. 

Within Quantum Monte Carlo, the Diffusion Monte Carlo is the method of choice for the 
calculation of ground-state properties of appropriate reference homogeneous systems, (the 
path integral method^ can be used to extend the study to non-zero temperatures degener- 
ate systems^), the most relevant example being the correlation energy of the electron gas 
obtained by Ceperley and Alder back in 1980^. This is even more so in the present days, 
since better wave-functions and optimization methods have been developed, better schemes 
to minimize finite-size effect have been devised, and vastly improved computational facilities 
are available. 

Recently, Zecca et ai.- have provided a Local Density functional for short-range pair 
potentials v{r) = eTic{fir)/r, whereas Paziani et ai.— have developed a Local Spin Density 
functional for the softened-core, long range case, v{r) = ei{{^r)/r. 

It is the purpose of this work to build on previous work-ii^ and provide the Radial 
Distribution Function (RDF), most notably the on-top value, i.e. its value at contact, at 
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given in Eq. ([2]). A brief account 



a zero radial distance, for the pair potential of Ref. 
of the results from this study has been presented in a recent communicationii. Aim of the 
present work is to give a complete and detailed account of the calculations that has been 
carried on for such a study. 

We performed fixed-nodes Diffusion Monte Carlo simulationsi^ where we used mod- 
ern techniques^ to optimize Slater- Jastrow wave-functions with backflow and three-body 



correlations^^ and Hellmann and Feynman (HFM) measures^ to calculate the RDF, par- 
ticularly the on-top value, which suffers from poor statistical sampling in its conven- 
tional histogram implementation. Twist-averaged boundary conditions'^ and RPA-based 
corrections'^ to minimize finite-size effects were not found essential for the RDF calculation. 

For the fully polarized and unpolarized fluid, we explored a range of densities and of the 
parameter /i. This required simulating several different systems. We also needed to evaluate 
and extrapolate out, for representative cases, time-step errors, population control bias, and 
size effects. We plan to explore intermediate polarizations in a future work. 

In the study we use two kinds of Jastrow-correlation-factors, one better for the near- 
Jellium systems and one better for the near-ideal systems. 

An important component of a computer experiment of a system of many particles, a 
fluid, is the determination of suitable estimators to measure, through a statistical average, a 
given physical quantity, an observable. Whereas the average from different estimators must 
give the same result, the variance, the square of the statistical error, can be different for 
different estimators. We compare the measure of the histogram estimator for the RDF with 
a particular HFM one. 

In ground state Monte Carlo simulationsi^ii^, unlike classical Monte Carlo simulation822""i2^ 
and path integral Monte Carlo simulations^, one has to resort to the use of a trial wave 
function^, \E'. While this is not a source of error, bias, in diffusion Monte Carlo simulation^ 
of a system of Bosons, it is for a system of Fermions, due to the sign problent^^. Another 
source of bias inevitably present in all three experiments is the finite size error. 

In a ground state Monte Carlo simulation, the energy has the zero-variance principle^^: 
as the trial wave function approaches the exact ground state, the statistical error vanishes. 
In a diffusion Monte Carlo simulation of a system of Bosons the local energy of the trial 
wave function, El(R.) = [i/^(R)]/\I'(R), where R denotes a configuration of the system of 
particles and H is the Hamiltonian, which we will here assume to be real, is an unbiased 
estimator for the ground state. For Fermions the ground state energy measurement is biased 
by the sign problem. For observables O which do not commute with the Hamiltonian the 
local estimator Ol(R) = [0^(R)]/^(R), is inevitably biased by the choice of the trial wave 
function. A way to remedy to this bias can be the use of the forward walking method^^i^ or 
the reptation quantum Monte Carlo method^l, to reach pure estimates. Otherwise this bias 
can be made of leading order 6'^ with 6 = (pQ — "^, where 0o is the ground state wave function, 



— ext 

introducing the extrapolated measure O = 2(0^)/ — {Ol) f^^^, where the first statistical 
average, the mixed measure, is over the diffusion Monte Carlo (DMC) stationary proba- 
bility distribution / and the second, the variational measure, over the variational Monte 
Carlo (VMC) probability distribution /vmc? which can also be obtained as the stationary 
probability distribution of a DMC without branchings^. 

One may follow different routes to determine estimators as the direct microscopic one, 
the virial route through the use of the virial theorem, or the thermodynamic route through 
the use of thermodynamic identities. This aspect of finding out different ways of calculating 
quantum properties in some ways resembles experimental physics. The theoretical concept 
may be perfectly well defined but it is up to the ingenuity of the experimentalist to find the 
best way of doing the measurement. Even what is meant by "best" is subject to debate. 
In an unbiased experiment the different routes to the same observable must give the same 
average. 

In this work we propose to use the Hellmann and Feynman theorem as a direct route for 
the determination of estimators in a diffusion Monte Carlo simulation. Some attempts in 
this direction have been tried before^^i^. The novelty of our approach is a different definition 
of the correction to the variational measure, necessary in the diffusion experiment, respect 
to Ref.— and the fact that the bias stemming from the sign problem does not exhaust all 
the bias due to the choice of the trial wave function, respect to Ref.—. 

The work is organized as follows: in Sec. [Til we introduce the fluid model; in Sec. IIIII 
we describe the Ewald sums technique to treat the long range pair-potential; in Sec. IIVI 
we describe the flxed-nodes Diffusion Monte Carlo (DMC) method; in Sec. |V]we describe 
several different ways to evaluate expectation values in a DMC calculation; in Sec. I VI I we 
describe the choice of the trial wave-function; in Sec. IVIII we deflne the RDF and describe 
some of its exact properties; the numerical results for the RDF are presented in Section 
IVIIIj Sec. IIXI is for flnal remarks. 

II. THE MODEL 

The Jellium is an assembly of A^ electrons of charge e moving in a neutralizing background. 
The average particle number density is n = N/Vt, where VL is the volume of the fluid. In the 
volume VL there is a uniform neutralizing background with a charge density pf, = —en. So 



that the total charge of the system is zero. 

In this paper lengths will be given in units of a = (47m/3)~^/^. Energies will be given in 
Rydbergs h'^/{2mal), where m is the electron mass and Qq = h^/{me^) is the Bohr radius. 

In these units the Hamiltonian of Jellium is 

1 ^ 

H = --2Y.K + v{R), (3) 

v=^{^y.T-^+y.^'+^o] , (4) 



r. \ ^-^ r, — To ^-^ ] 



where R = (ri, r2, . . . , r^) with Fj the coordinate of the ith electron, r^ = a/ao, and vq a 
constant containing the self energy of the background. 

The kinetic energy scales as 1/r^ and the potential energy (particle-particle, particle- 
background, and background-background interaction) scales as l/vg, so for small r^ (high 
electronic densities), the kinetic energy dominates and the electrons behave like an ideal gas. 
In the limit of large r^, the potential energy dominates and the electrons crystallize into a 
Wigner crystal^i^. No liquid phase is realizable within this model as the pair-potential has 
no attracti ye p arts even though a superconducting state^^ may still be possible (see chapter 



8.9 of Ref. 



M). 



A. Modified long range pair-potential 

The fluid model studied in this work is obtained modifying the Jellium by replacing 
the 1/r Coulomb potential between the electrons with the following long range bare pair- 
potential^ 

erf(/xr) 
v^,{r) = —— , (5) 

whose Fourier transform is 

An _ k^ 
Mk) = p e-3^ . (6) 

When fi — )■ oo, we recover the standard Jellium model; in the opposite limit /i — ?■ 0, we 
recover the non-interacting electron gas. Notice that t'^ is a long range pair-potential with 
a penetrable core, v^{0) = 2^j ^. So /i controls the penetrability of two particles. For 
this kind of system it is lacking a detailed study of the RDF. In this work we will only be 
concerned about the fluid phase. 
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v{r) = / —— e-'^^'vik) , v{k) = / dre'^^^r) . (7) 



III. EWALD SUMS 

Periodic boundary conditions are necessary for extrapolating results of the finite system 
to the thermodynamic limit. Suppose the bare pair-potential, in infinite space, is v{r), 

dk 
(2^3 

The best pair-potential of the finite system is given by the image potential 

y/(r) = ^t;(|r + L|)-f;(0)/f]. (8) 

L 

where the L sum is over the Bravais lattice of the simulation cell L = (m^X, myL, rrizL) 
where mx,my,mz range over all positive and negative integers and Q = L^. We have also 
added a uniform background of the same density but opposite charge. Converting this to 
/c-space and using the Poisson sum formula we get 

vi{r) = ^i2^'{k)e-'''-\ (9) 

k 

where the prime indicates that we omit the k = term; it cancels out with the background. 
The k sum is over reciprocal lattice vectors of the simulation box kn = (27™^,/-^, 2'Kny/L, 
2Tinz/L) where nx,ny,nz range over all positive and negative integers. 

Because both sums, Eq. ([8]) and Eq. iQ, are so poorly convergent^i we follow the scheme 
put forward by Natoli et ai.— for approximating the image potential by a sum in /c-space 
and a sum in r-space, 

^,(r) = Y,Vs{\r + L\)+J2 ^liky'^'^ " m/^ , (10) 

L |k|<fcc 

where Vs{t) is chosen to vanish smoothly as r approaches Tc, where r^ is less than half of 
the distance across the simulation box in any direction. If either r^ or kc go to infinity 
then Va ^ vj. Natoli et ai. show that in order to minimize the error in the potential, it is 
appropriate to minimize x^ = Inl'^ii^) ~ '^'a{r)]'^ dr/Q. And choose for Vs{r) an expansion in 
a fixed number of radial functions. This same technique has also been applied to treat the 
Jastrow-correlation-factor described in section IVI A[ 

Now let us work with N particles of charge e in a periodic box and let us compute the 
total potential energy of the unit cell. Particles i and j are assumed to interact with a 



potential e^v{rij) = e^f (|rj — rj\). The potential energy for the A^ particle system is 

y = ^eV(r.,) + ^e2^M, (11) 

i<j i 

where vm = \ linir_s.o[t'/(r) — t'(r)] is the interaction of a particle with its own images; it is a 
Madelung constant^ for particle i interacting with the perfect lattice of the simulation cell. 
If this term were not present, particle i would only see A^ — 1 particles in the surrounding 
cells instead of N . 

IV. THE FIXED-NODES DIFFUSION MONTE CARLO (DMC) METHOD 

Consider the Schrodinger equation for the many-body wave-function, 0(R, t) (the wave- 
function can be assumed to be real, since both the real and imaginary parts of the wave- 
function separately satisfy the Schrodinger equation), in imaginary time, with a constant 
shift Et in the zero of the energy. This is a diffusion equation in a 3A^-dimensional spaced. 
If E-x is adjusted to be the ground-state energy, E'q) the asymptotic solution is a steady 
state solution, corresponding to the ground-state eigenfunction 0o(R) (provided 0(R, 0) is 
not orthogonal to 0o)- 

Solving this equation by a random- walk process with branching is inefficient, because 
the branching rate, which is proportional to the total potential V^(R), can diverge to -|-oo. 
This leads to large fluctuations in the weights of the diffusers and to slow convergence when 
calculating averages. However, the fluctuations, and hence the statistical uncertainties, can 
be greatly reduced^ by the technique of importance sampling^. 

One simply multiplies the Schrodinger equation by a known trial wave-function ^^(R) 
that approximate the unknown ground-state wave-function, and rewrites it in terms of a 
new probability distribution 

/(R,t)=</)(R,t)vl/(R), (12) 

whose normalization is given in Eq. ( JAll) . This leads to the following diffusion equation 

-AVV(R, t) + \El[B.) - i5;T]/(R, t) + AV ■ [/(R, t)F(R)] . (13) 



5/(R,^) 



dt 

Here A = f? j{2m\ t is the imaginary time measured in units of h, El(R) = [iir\l/(R)]/\E'(R) 
is the local energy of the trial wave-function, and 

F(R) =Vln^2^R) . (14) 

8 



The three terms on the right hand side of Eq. (fT3|) correspond, from left to right, to diffusion, 
branching, and drifting, respectively. 

At sufficiently long times the solution to Eq. ( TTS]) is 

fiR,t)^No^(R)M^)eM-iEo-ET)t] , (15) 

where A'o = / 0o(R)0(R, 0) rfR. If Et is adjusted to be Eq, the asymptotic solution is 
a stationary solution and the average (£'i(R))/ of the local energy over the stationary 
distribution gives the ground-state energy Eq. If we set the branching to zero El(II) = Et 
then this average would be equal to the expectation value J ^(R)i/^(R) rfR, since the 
stationary solution to Eq. ( TT3l) would then be / = /vmc = ^I/^. In other words, without 
branching we would obtain the variational energy of \1/, rather than Eq, as in a Variational 
Monte Carlo (VMC) calculation. 

The time evolution of /(R, t) is given by 

/(R', t + r) = j rfRG(R', R; r)/(R, t) , (16) 

where the Green's function G(R', R; r) = *(R')(R'| exp[-r(/7 - Er)]|R)*~HR) is a tran- 
sition probability for moving the set of coordinates from R to R in a time r. Thus G 
is a solution of the same differential equation, Eq. (fT3l) . but with the initial condition 
^(R , R; 0) = (5(R — R). For short times r an approximate solution for G is 

G(R', R; r) = (4vrAr)"=^^/2e^^'-^-^-^WI'/^^-e"-{[^-(^)+^-(^')l/2-^-> + Oir') . (17) 

To compute the ground-state energy and other expectation values, the A^-particle distribu- 
tion function /(R, t) is represented, in diffusion Monte Carlo, by an average over a time 
series of generations of walkers each of which consists of a fixed number of n^ walkers. A 
walker is a pair (Rq,,cJq,), a = l,2,...,n^, with Rq, a 3A^-dimensional particle configu- 
ration with statistical weight Ua- At time t, the walkers represent a random realization 
of the A^-particle distribution, /(R, t) = J2a=i'^a^(^ " ^a)- The ensemble is initial- 
ized with a VMC sample from /(R, 0) = \1/^(R), with w° = 1/ra^ for all a. Note that 
if the trial wave-function were the exact ground-state then there would be no branching 
and it would be sufficient n^ = 1. A given walker (R*,a;*) is advanced in time (dif- 
fusion and drift) as R*"*""^ = R* -|- x + ArVln\&^(R*) where x is a normally distributed 
random 3A^-dimensional vector with variance 2Ar and zero mean^S. In order to satisfy de- 
tailed balance we accept the move with a probability v4(R, R'; r) = min[l, H^(R, R')]) where 



W{R,R') = [G{R,R';t)<I/\R')]/[G{R',R;t)-^^{R)]. This step would be unnecessary if G 
were the exact Green's function, since W would be unity. Finally, the weight cu^ is replaced 
by ui+^ = uji^ui (branching), with Aw^ = exp{-r[(Ei(R^^) + El{R'+-))/2 - Et]}. 

However, for the diffusion interpretation to be valid, / must always be positive, since it 
is a probability distribution. But we know that the many-fermions wave- function 0(R, t), 
being antisymmetric under exchange of a pair of particles of the parallel spins, must have 
nodes, i.e. points R where it vanishes. In the fixed- nodes approximation one restricts the 
diffusion process to walkers that do not change the sign of the trial wave-function. One can 
easily demonstrate that the resulting energy, {El{R)) f, will be an upper bound to the exact 
ground-state energy; the best possible upper bound with the given boundary condition^^. 

A detailed description of the algorithm used for the DMC calculation can be found in 
Ref. 28. 



V. EXPECTATION VALUES IN DMC 

In a DMC calculation there are various different possibilities to measure the expectation 
value of a physical observable, as for example the RDF. If (O)/ is the measure and (. . .) f 
the statistical average over the probability distribution / we will, in the following, use the 
word estimator to indicate the function O itself, unlike the more common use of the word to 
indicate the usual Monte Carlo estimator X]i=i ^il-^ of ^h^ average, where {Oi] is the set 
obtained evaluating O over a finite number M of points distributed according to /. Whereas 
the average from different estimators must give the same result, the variance, the square of 
the statistical error, can be different for different estimators. 

1. The local estimator and the extrapolated measure 

To obtain ground-state expectation values of quantities O that do not commute with 
the Hamiltonian we introduce the local estimator Ol(R) = [0\E'(R)]/\E'(R) and then com- 
pute the average over the DMC walk, the so called mixed measure, O = (Ol(R))/ = 
/ (/)o(R)0\E'(R) dR/ J (/)o(R)\&(R) dR. This is inevitably biased by the choice of the trial 
wave-function. A way to remedy to this bias is the use of the forward walking method^^^^^ or 
the reptation quantum Monte Carlo method^i to reach pure estimates. Otherwise this bias 
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can be made of leading order 6'^ , with 5 



\1/, introducing the extrapolated measure 



■p^TCxt 



20 - O 



(18) 



where O 



{Ol) f^mc is the variational measure. If the mixed measure equals the varia- 
tional measure then the trial wave-function has maximum overlap with the ground-state. 



2. The Hellmann and Feynman measure 

Toulouse et ai.— 12^ observed that the zero-variance property of the energy^i can be ex- 
tended to an arbitrary observable, O, by expressing it as an energy derivative through the 
use of the Hellmann- Feynman theorem. 

In a DMC calculation the Hellmann-Feynman theorem takes a form different from the 
one in a VMC calculation. Namely we start with the eigenvalue expression {H^ — E^)'$^ = 
for the ground-state of the perturbed Hamiltonian i/^ = H + XO, take the derivative with 
respect to A, multiply on the right by the ground-state at A = 0, (po, and integrate over the 
particle coordinates to get 



rfR0o(i^ -^ 



dRi 



^' 



(19) 



d\ J ""V d\ d\ 
Then we notice that due to the Hermiticity of the Hamiltonian, at A = the left hand side 
vanishes, so that we geir^ 



JdR(f)oO^^ 



/t/R0o*^ 



OE^ 



A=0 



dX 



(20) 



A=0 



This relation holds only in the A — ;■ limit unlike the more common fornti^ which holds for 



30 



any A. Also it resembles Eq. (3) of Ref. 

Given E^ = / dR(/)o(R)i/^^^(R)/ /dR0o(R)^^(R) the "Hellmann and Feynman" 
(HFM) measure in a DMC calculation is 

dE^ 

A=0 



O 



■HFM 



d\ 



1^1 



(0^(R)); + (A02(R))^ + (A0^(R)) 



/ 



(21) 



The a correction is^ 



A02(R) 






ElIR) 






(22) 



This expression coincides with Eq. (18) of Ref. Il5|. In a VMC calculation this term, usually, 
does not contribute to the average, with respect to /vmc = ^^5 due to the Hermiticity of 
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the Hamiltonian. This is of course not true in a DMC calculation. We will then define a 
Hellmann and Feynman variational (HFMv) estimator as 0^^^^"" = Ol(R) + A0£(R). The 
/5 correction is^i 

AOliR) = [E,{R) - Eo]^ , (23) 
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by a factor of one half. This term 



where Eq = E''^^^. Which differs from Eq. (19) of Ref. 

is necessary in a DMC calculation not to bias the measure. The extrapolated Hellmann and 

Feynman measure will then be 



^rHFM-ext _^tHFM , ^HFMv 



O '-^^=20' -(0^^^^)/_. (24) 

Both corrections a and /3 to the local estimator depends on the auxiliary function, \E'' = 
d'^^/dX\\=o. Of course if we had chosen \[f^=°, on the left hand side of Eq. ( 1211) . as the exact 
ground state wave-function, (pQ, instead of the trial wave- function, then both corrections 
would have vanished. When the trial wave-function is sufficiently close to the exact ground 
state function a good approximation to the auxiliary function can be obtained from first 
order perturbation theory for A ^ 1. So the Hellmann and Feynman measure is affected by 
the new source of bias due to the choice of the auxiliary function independent from the bias 
due to the choice of the trial wave-function. 

It is convenient to rewrite Eqs. (p2l) and (123|) in terms of the logarithmic derivative 
Q(R) = ^'(R)/*(R) as follows 

1 ^ 
A02iR) = -- J2[VlQ{R) + 2v,(R) ■ Vr,g(R)] , (25) 

''^'^ k=i 

AOf(R) = [Ei(R)-E]g(R), (26) 

where v/c(R) = Vrj. ln\E'(R) is the drift velocity of the trial wave- function. For each observ- 
able a specific form of Q has to be chosen. 

VI. TRIAL WAVE-FUNCTION 

We chose the trial wave-function of the Bijl-Dingle-Jastrow^ or product form 

^(R) oc D(R) exp [ - ^ u{rij) j . (27) 
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The function -D(R) is the exact wave- function of the non-interacting fermions (the Slater 
determinant) and serves to give the trial wave-function the desired antisymmetry 

DiR) = -^ det(^+„)-^ det(^-^) , (28) 

where for the fluid phase (Pnm = c*''"'^'"5o-,„,o-/V" with k„ a reciprocal lattice vector of 
the simulation box such that |k.„| < kp, a the ^-component of the spin (±1/2), r^ the 
coordinates of particle m, and am its spin z-component. For the unpolarized fluid there are 
two separate determinants for the spin-up and the spin-down states because the Hamiltonian 
is spin independent. For the polarized fluid there is a single determinant. For the general 
case of A'^.i. spin- up particles the polarization will be (^ = (A^+ — NJ)/N and the Fermi wave- 
vector for the spin-up (spin-down) particles will be kp = {l^CY^^kp with kp = (37r^n)^'^ = 
(97r/4)^/'^/(aors) the Fermi wave- vector of the paramagnetic fluid. On the computer we fill 
closed shells so that N^j is always odd. We only store k„ for each pair (k,„, — k^) and use 
sines and cosines instead of exp(zk„ ■ Fj) and exp(— zk„ ■ Tj). 

The second factor (the Jastrow factor) includes in an approximate way the effects of 
particle correlations, through the "Jastrow-correlation-factor" , u{r), which is repulsive. 

A. The Jastrow-correlation-factor 

Neglecting the cross term between the Jastrow and the Slater determinant in Eq. f lA6p 
(third term) and the Madelung constant, the variational energy per particle can be approx- 
imated as follows, 

k 

I ' 

^^^Ak-kM(/i;)M(/i;')(pk+k'P-kP-k')/ + --. , (29) 

k,k' 

where ep = (3/5)A ^^^ Act(A;^)^/A is the non-interacting fermions energy per particle, u{k) 
is the Fourier transform of the Jastrow-correlation-factor M(r), v^{k) = Air expl—k"^ / 4fi^) / k'^ 
is the Fourier transform of the bare pair-potential, S{k) is the static structure factor for 
a given u(r) (see Sec. IVll 31) . pk = X]j=i6xp(«k ■ r^) is the Fourier transform of the total 
number density p(r) = X]j^('^ ~ ^i)^ ^^^ ^^^ trailing dots stand for the additional terms 
coming from the exclusion of the j = k term in the last term of Eq. (IA6p . Next we make 
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the Random Phase Approximation:^ and we keep only the terms with k + k' = in the last 
term. This gives 



+ — ^Ue^v^ik) - 2\eu{k)][S{k) - 1] - 2n\[ku{k)]'^S{k)\ + ... . (30) 

k 

In the limit fc — )> we have to cancel the Coulomb singularity and we get u'^{k) = 



k 



me^Vfj, 



{k) / [h^nk'^) ~ [{Aixe^ / k"^) / {hujp)]^ (where ujp = -^f Aime^ / m is the plasmon frequency) 
or in adimensional units 

u{k) = W^— , small k . (31) 

This determines the correct behavior of u{k) as A; — )■ or the long range behavior of u{r) 

u{r) = J^- , large r . (32) 

Now to construct the approximate Jastrow-correlation-factor, we start from the expres- 
sion 

e = eF + ^J2le%{k)-AXeu{k)][S{k)-l] , (33) 



k 

\42 



and use the following perturbation approximation, for how S{k) depends on u{k)^, 

1 1 



S(k) S^(fc) + '^"«'*) • (^^' 

where A and B are constant to be determined and S^{k) the structure factor for the non- 
interacting fermions (see Eq. (!62|) ). which is S^ = ^^ S^^ with 

Jj''-'-'"^' (35) 

— else 

n 

where no- = N^/Vt and i/o- = k/{2kp). 

Minimizing e with respect to u{k), we obtain^ 



Bnu{k) 



S'^ik) 



Bne^v^{k) 



1/2 



(36) 



_S^{k) XAk^ 

This form is optimal at both long and short distances but not necessarily in between. In 
particular, for any value of C, the small k behavior oiu{k) is y2r7/3^AS(47r/A;^) which means 
that 



2r 1 
«W = \/^-' larger. (37) 
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The large k behavior oiu{k) is {ts/ A)v^{k)/k'^^ for any value of C, which in r space translates 
into 



du{r) 



dr 



— 7 U ^- oo 
2A ^ 



r=0 







(38) 



fi finite 

In order to satisfy the cusp condition for particles of antiparallel spins (any reasonable 
Jastrow-correlation-factor has to obey to the cusp conditions (see Ref. Il3| Section IVF) 
which prevent the local energy from diverging whenever any two electrons {fi = oo) come 
together) we need to choose A = 1, then the correct behavior at large r f lHT]) is obtained 
fixing B = 2 (see Note 1441) . We will call this Jastrow J7i in the following. 

It turns out that, at small /x, but not for the Coulomb case, a better choice is given by^ 

1/2 



2nu{k) 



S^{k) 



S^{k) 



2ne Vfj,{k) 



(39) 



which still has the correct long ( !37l) and short (|38l) range behaviors. We will call this Jastrow 
^2 in the following. This is expected since, differently from J7i, 1/2 satisfies the additional 
exact requirement lim^_i.o u{r) = 0, as immediately follows from the definition (l39l) . Then, as 
confirmed by our results (see Sec. IVIIIEp ). at small /i (and any r^), the trial wave-function 
is expected to be very close to the stationary solution of the diffusion problem. 



B. The backfiow and three-body correlations 



As shown in Appendix 13 the trial wave-function of Eq. ( 1271) can be further improved by 
adding three-body (3B) and backfiow (BF) correlations^"^ as follows 

N 



^(R) = Z)(R)exp 



5^^(r,,)-5^G(/)-G(/) 



i<j 



1=1 



Here 



D{R) 



V^ 



det{(p 



n,ni/ 



v/AM 



det(<^„. 



with 7?° 



^2K7T,-Xj7 



(^o-jn.cr/v f^ and x^ quasi-particle coordinates defined as 



TV 



Xi = rj + ^ r]{rij){ri - r^] 



(40) 



(41) 



(42) 
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The displacement of the quasi-particle coordinates Xj from the real coordinate r j incorporates 
effects of hydrodynamic backfiow^, and changes the nodes of the trial wave- function. The 
backflow correlation function T][r), is parametrized as^ 

Vir) = \b l+'^l , , (43) 

which has the long-range behavior ~ 1/r^. 

Three-body correlations are included through the vector functions 

N 

G{^) = J2^ir^J){r^-r,). (44) 

We call ^(r) the three-body correlation function which is parametrized as^ 

e(r) = aexp{-[(r-6)c]2} . (45) 

To cancel the two-body term arising from G(/) ■ G(/), we use u{r) = u{r) — 2^^(r)r^ 

The backflow and three-body correlation functions are then chosen to decay to zero with 
a zero first derivative at the edge of the simulation box. 

C. Optimized parameters 



Optimizing the trial wave-function (see Ref. Il3| Section VII) is extremely important for 
a fixed-nodes DMC calculation as, even if the Jastrow-correlation-factor is parameter free, 
the backflow changes the nodes. We carefully studied how the RDF depends on the quality 
of the trial wave-function choosing a simple Slater determinant (S) (Eq. f l27|) without the 
Jastrow factor), a Slater- Jastrow (SJ) (Eq. (1271) ). and a Slater- Jastrow with the backflow 
and three-body corrections (SJ+BF-F3B) (Eq. (I40|l ). 

In Table [T] we report the optimized parameters for the backflow and three-body correlation 
functions for a system of A^ = 54 and C = at various r^ and /i. We have used these values 
of the parameters in all subsequent calculations, unrespective of the value of C,. 

In Fig. [T]we show the optimized rj and ^ for A^ = 54, C = 0, r^ = 10. The optimization of 
the 7 parameter dependent trial wave-function gives a backflow correlation rj ordered in /x 
but a three-body correlation ^ erratic in yU. As one moves away from the Coulomb yU — )■ oo 
case the system of particles becomes less interacting and the relevance of the backflow 
and three-body correlations diminishes. This is supported by the fact that at /x = 4, 2, 1, 
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TABLE I. Optimized variational parameters of backflow and three-body correlation functions for 
A^ = 54 and C = ^-iid various combinations af r^ and fi. 



Tb 


M 


As 


SB 


r-B 


wb 


a 


b 


c 


10 


1/2 


- 


- 


- 


- 


- 


- 


- 


10 


1 


8.408d-4 


1.658d+2 


-1.383d-3 


3.168 


0.447 


-0.212 


1.036 


10 


2 


7.189d-5 


9.793d+2 


9.478d-6 


0.446 


1.379d-fl 


-3.688 


0.450 


10 


4 


1.116d-4 


6.522d+2 


-2.553d-5 


0.179 


5.981d+l 


-4.773 


0.462 


10 


oo 


0.781 


-0.499 


0.324 


2.958 


0.514 


0.327 


1.358 


5 


1/2 


- 


- 


- 


- 


- 


- 


- 


5 


1 


- 


- 


- 


- 


- 


- 


- 


5 


2 


2.768d-2 


-0.420 


0.893 


-0.673 


1.322d-|-6 


-9.003 


0.408 


5 


4 


0.331 


-0.680 


1.467 


1.442 


2.729d+l 


-2.607 


0.659 


5 


oo 


0.161 


-0.585 


0.335 


0.841 


0.802 


-7.310d-2 


1.344 


2 


1/2 


- 


- 


- 


- 


- 


- 


- 


2 


1 


- 


- 


- 


- 


- 


- 


- 


2 


2 


- 


- 


- 


- 


- 


- 


- 


2 


4 


5.272d-2 


-1.616 


1.732 


1.687d-2 


804.135 


-2.875 


0.847 


2 


oo 


5.018d-2 


-1.221 


0.393 


0.681 


1.655 


-0.596 


1.229 




1/2 


- 


- 


- 


- 


- 


- 


- 




1 


- 


- 


- 


- 


- 


- 


- 




2 


- 


- 


- 


- 


- 


- 


- 




4 


1.187d-2 


-6.834 


0.495 


1.295 


0.186 


0.489 


4.739 




oo 


2.1945d-2 


-3.086 


0.320 


1.631 


0.306 


0.367 


2.467 



in correspondnce of the erratic behavior, the effect of the three-body correlations on the 
expectation value of the energy is irrelevant. 
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FIG. 1. Shows the optimized correlation functions r] and ^ for N = 54, C = 0, and r^ = 10 and 
different values of /i. 

VII. THE RADIAL DISTRIBUTION FUNCTION (RDF) 

The main purpose of the present work is to determine the radial distribution function 
(RDF) of our fluid model through the DMC calculation. 



1. Definition of the radial distribution function 

The spin-resolved RDF is defined as^*^ 
5'<7,<7'(r,r ) = 



ra„ r m„/ r 



''(T \^ I '^CF' 



N 



rirrir 



E^<^'-'^( 



r - Ti 



(46) 
(47) 



. i=l 



where here, and in the following, (...) will denote the expectation value respect to the 
ground-state. Two exact conditions follow immediately from the definition: i. the zero- 
moment sum rule 



^ / drdr' n„{r)n^{r')[g„y{r,r') - 1] = -N 



(48) 
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also known as the charge (monopole) sum rule in the sequence of multipolar sum rules in 
the framework of charged fluids^, ii. gcr,a{T^, r) = due to the Pauli exclusion principle. 

For the homogeneous and isotropic fluid no-(r) = Nf^/Vt where No- is the number of 
particles of spin a and g^^a' depends only on the distance r = |r — r |, so that 

9-Ar) = 4^^2jv^^^, ( H KoA',o,5{r - n,) \ . (49) 

The total (spin-summed) radial distribution function will be 

n^ z — / 

^ ^ ^^ 9^Ar) + f^l 9-Ar) + ^^+,-(r) ■ (50) 



2. From the structure to the thermodynamics 

As it is well known the knowledge of the RDF gives access to the thermodynamic prop- 
erties of the system. The mean potential energy per particle can be directly obtained from 
g{r) and the bare pair-potential v^{r) as follows 



6p 



E^/^r^VW[^.,.'W-i], (51) 



where we have explicitly taken into account of the background contribution. Suppose that 
ep{rs) is known as a function of the coupling strength r^. The virial theorem for a system 
with Coulomb interactions (foo('") = I/'") gives N{2ek + Cp) = 3PQ with P = —d{NeQ)/dQ 
the pressure and Cq = 6^ + Cp the mean total ground-state energy per particle. We then find 

e,(r.) = 2eo(r.) + r.^ = ^^Heo(r.)] , (52) 



s "" s 



which integrates to 

1 T" 
eo(r-.) = eF + — / dr'^r'^ep{r'^) . (53) 

n Jo 

We can rewrite the ground-state energy per particle of the ideal Fermi gas, in reduced 

units, as 

fdnA^ 3 , , 1 



s 
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where 4>n{C) = (1 ~ Q'^/^ + (1 + C)""^^- And for the exchange potential energy per particle in 
the Coulomb case 

which follows from Eq. ( 15T1) and Eqs. (l59l) -( !6Ql) . The expression for finite n can be found in 



Ref. 



id (see their Eqs. (15)-(16)). 



3. Definition of the static structure factor 

If we introduce the microscopic spin dependent number density 

N 



j=l 



and its Fourier transform Pk,o-; then the spin-resolved static structure factors are defined as 
S'cr,(j'(k) = {pk,aP-'k,a')/N, which, for the homogeneous and isotropic fluid, can be rewritten 
as 

S.Ak) = -S.y + ^^^^ [[g.y{r) - l]e-^^-^dr + ^^^^(27r)3<5(k) , (57) 

n n J n 

From now on we will ignore the delta function at k = 0. The total (spin-summed) static 
structure factor is S = Ylcra' ^f^,^'- ^^^ ^'^ ^^^ charge sum rule f HSjl we must have 
linifc^o S{k) = 0. In Sec. IVII B 21 we will show that the small k behavior of S{k) has 
to start from the term of order k"^. 



A. Analytic expressions for the non-interacting fermions 

Usually Qa^cr' is conventionally divided into the (known) exchange and the (unknown) 
correlation terms 

9a,a' = 9a,a' + 9a,a' > (58) 

where the exchange term corresponds to the uniform system of non-interacting fermions. 
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1. Radial distribution function 

We thus have (from the definition of the RDF (146|) and using Slater determinants for the 
wave-function) 

9lA^) = 1 , (59) 



.X r^\_l_ 



9a A^) 



3ji(/c^r) ^^ 
kpT 



(60) 



where ji{x) = [sin(x) — xcos(x)]/x^ is the spherical Bessel function of the first kind and 
(kp)^ = Gvr^rio- is the Fermi wave-number for particles of spin a. 



2. Static structure factor 

Again we will have the splitting S^y = S^^, + S^^, into the exchange and the correlation 
parts. So that for the non-interacting fermions we get 



SlJk) = 0, (61) 

Si,W.^-!ie,2.J-.)^(l-J,)^(2 



Zkp 



rifj \ ^ /c > 2kp 

n I 3_fc 1_ ( _k_\ U ^ OLo- 



4 fcf, 16 \^ k'p 



(62) 



where 9(x) is the Heaviside step function. 



B. RDF sum rules 

Both the behavior of the RDF at small r and at large r has to satisfy to general exact 
relations or sum rules. 
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1. Cusp conditions 

When two electrons (/i = oo) get closer and closer together, the behavior of gayif') is 
governed by the exact cusp conditions^ 



d 

^s9.Ar) 

d 
ar 



= 0, 




r->0 


3d' , , 




r->0 ^0 





(63) 
(64) 
(65) 



where in the adimensional units Qq — )■ l/r^. For finite /i we only have the condition 5'cr,<T(0) = 
due to Pauli exclusion principle. 

2. The Random Phase Approximation (RPA) and the long range behavior 
of the RDF 

Within the linear density response theory^i^ one introduces the space-time Fourier 
transform, x(k, cj), of the linear density response function. Which is related through 
the fiuctuation dissipation theorem, S'(k, w) = — (2^/r2)9(ci;)Imx(k, w), to the space-time 
Fourier transform, S'(k, uo) (dynamic structure factor), of the van Hove correlation function^, 
{p{r,t)p{O,0))/n, where p(r,t) = exp{iHt/h)p{r)exp{—iHt/h). 

In the Random Phase Approximation (RPA) we have^ 

71 — V = ~Ti V ~ e'v^{k) , (66) 

where Xo is the response function of the non-interacting Fermions (ideal Fermi gas), known 
as the Lindhard susceptibility^. This corresponds to taking the "proper polarizability" (the 
response to the Hartree potential) equal to the response of the ideal Fermi gas^. The RPA 
static structure factor is then recovered from the fiuctuation dissipation theorem as follows 

SRPAik) = / — lmxRPAik,uj) . (67) 

n Jo TT 



where 



ImxRPA = jz ^T-5 ^2 ° / 2- T ^^ ' (^^) 
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The small k behavior of the RPA structure factor is exact^. One finds 

hk'^ 

SRPAik) = , k<^kF , (69) 

2mUp 

where Up = A/47me^/m is the plasmon frequency—. This is also known as the second-moment 
sum rule for the exact RDF and can be rewritten as n j dv r'^[g{r) — 1] = '-6{h/2mUp). We 
can then say that g{r) — l has to decay faster than r~^ at large r. The fourth- moment (or com- 
pressibility) sum rule links the thermodynamic compressibility, x = [nd{n'^deo/dn)/dn\~^,— 
to the fourth-moment of the RDF. For the equivalent classical system it is well known that 
the correlation functions have to decay faster than any inverse power of the distance^>^>^ 
(in accord with the Debye-Hiikel theory). We are not aware of the existence of a similar 
result for the zero temperature quantum case. 

VIII. RESULTS OF THE CALCULATION 

We considered fourty systems corresponding to r^ = 1,2,5,10, /i = oo,4, 2, 1, 1/2, ( = 
0, 1. For each system we calculated the RDF using the histogram estimator in a variational, 
mixed, and extrapolated measure and a particular HFM measure. Before starting with the 
simulations we determined the optimal values for the time step r and the number of walkers 
Uw for each density. 

A. Extrapolations 

For the Coulomb case, yU — )■ oo, we made extrapolations in time step r and number of 
walkers n^ for each value of r^ within our DMC simulations. Given a relative precision 
6eo = Aeo/Cp, where eo = {El) f/N , Acq is the statistical error on cq, and e^ is the exchange 
energy per particle (see Eq. (15^ ). we set as our target relative precision 5^^ = 10~^%. 

1. In time step 

Our results are summarized in Table HIl As the characteristic dimension of one particle 
diffusing walk is cr = v2Ar or -J^^rjrl in adimensional units, this has to remain of the order 
of the mean nearest neighbor separation a which is chosen to be a constant in our units. 
Then we expect that at lower r^ one needs to choose smaller time steps r. For this reason 
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TABLE II. Extrapolation in time step for A^ = 66 unpolarized electrons (/U = oo) at a fixed number 
of riuj = 600 walkers with a trial wave-function of the S J type. We run the simulation for 3 different 
time steps and did a linear fit of the (r, cq) data, cq = a + br. The optimal r is the largest one 
compatible with the target precision. 

rs a b x^ optimal r 

10 -0.107456(7) 0.00010(2) 0.9 0.09 

5 -0.153352(4) 0.00024(3) 0.1 0.07 

2 -0.00416(8) 0.003(2) 4.4 0.01 

1 1.14579(7) 0.032(9) 1.1 0.003 

TABLE III. Extrapolation in number of walkers for A = 66 unpolarized electrons (/i = oo) with 
a time step r = 0.1 for r^ = 10,5, r = 0.05 for r^ = 2, and r = 0.01 for r^ = 1 with a trial 
wave-function of the SJ type. We run the simulation for 4 different numbers of walkers and did 
a linear fit of the (1/n^, eo) data, sq = a + b/riuj. The optimal n^ is the smallest one compatible 
with the target precision. 



10 -0.107443(3) 

5 -0.153329(6) 

2 -0.004036(6) 

1 1.14609(6) 



we chose different time steps in the simulations of the Table: r = 0.5, 0.1, 0.05 for r^ = 10, 
T = 0.3,0.1,0.05 for r, = 5, r = 0.05,0.03,0.005 for r, = 2, and r = 0.01,0.005,0.001 for 
Tj, = 1. Note that, at fixed Vg, the statistical errors increase as the time step diminishes. 



b 


x' 


optimal n, 


0.0032(4) 


0.1 


354 


0.0044(7) 


0.2 


243 


0.0026(7) 


0.2 


56 


0.01(1) 


1.2 


40 



2. In the number of walkers 

Our results are summarized in Table IIIII The fluctuations of the statistical weight of 
a walker depend on the fluctuations of the local energy, i.e. by the quality of the trial 
wave-function. The quality of the trial wave-function worsens as Vg becomes larger (for the 
strongly correlated system), and one expects that the necessary number of walkers increases. 
This is in agreement with the results of the Table. Note that, at flxed r^, the statistical 
errors increase as the number of walkers diminishes. 
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B. Effect of backflow and three-body correlations 

In Fig. [2] we sliow tlie mixed measure of tlie RDF calculated in DMC for A^ = 162, ( = 
0,/i = oOjTj, = 10 with different kinds of trial wave- functions. Of course in a VMC calcula- 
tion using the Slater determinant wave-function gives us g^^,, the RDF of the ideal gas (see 
Eqs. dSnD-dSO])). 
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FIG. 2. Shows the mixed measure of the RDF calculated in DMC for N = 162, ( = 0, fi = oo, r^ 
10 with a S, SJ, SJ+BF+3B trial wave- function. 



In Fig. in] we show the difference between the RDF calculated with the SJ wave-function 
and the one calculated with the SJ-I-BF-I-3B wave-function using the variational, the mixed, 
and the extrapolated measure. 
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FIG. 3. Shows the difference between the RDF calculated with the SJ wave-function and the one 
calculated with the SJ+BF+3B wave-function using the variational, the mixed, and the extrapo- 
lated measure. The results are for N = 162, C = 0,;U = oo. On the left the like RDF is used at 
Tg = 10, on the right the unlike RDF at r^ = 1. 



With the extrapolated measure the results from the SJ computation differs by less than 
0.005 from the ones from the SJ+BF+3B. We then decided to perform our subsequent 
calculations using the SJ trial wave-function. 



C. Size effects 

In order to estimate the size effects on the RDF calculation we performed a series of VMC 
calculation with the SJ wave-function on an unpolarized system with different number of 
particles. The results (see Fig. H]) show that the size dependence mainly affects the long 
range behavior of the RDF and the on-top value for the unlike one. 
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FIG. 4. Shows the difference between the RDF of two systems of electrons (/x = oo) at r^ = 10 and 
(^ = with different sizes Ni and N2. The RDF are calculated in VMC with the SJ wave- function. 
On the left the difference of the like RDF is shown and on the right the difference of the unlike 
RDF is shown. 

In the simulation the RDF is defined on r G [0, rmax] with Vmax = L/2 where L = Q^^^ = 
{AttN/SY'^ is the size of the simulation box. To minimize size effects we chose to perform 
our RDF calculation with A^ = 162 in the unpolarized case and A^ = 147 in the polarized 
case. 

D. The HFM measure 



From the definition (1491) . we can write the RDF as 

(J,,,.(r,R)) 



9a,a'[r) 



Qn^jTla' 



Since the operator I^ry is diagonal in coordinate representation then I^ry 
eating with Qr the solid angle spanned by the r vector, we can write 






(70) 



{Ia,a')L- Indi- 



(71) 
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which is the usual histogram estimator—. Following Toulouse^^ we choose for Q the following 
expression 

Q'Ar,B) = -f^Y.^„M„^y^^, (72) 

so that (using the identities Vr^^, 1/1 r- Fiji = -47r5(r-r.y) and Vr,/(rfcj) =VrkJijkj)[5ik- 
(5jj], for a given function /) the first term in Eq. fl25|l exactly cancels the histogram estimator 
Ir, r,i ■ Then the HFMv estimator is 



'-rj.cT 



1 v^ . . ._. r 



-I^ 2^ 5.,.,5.',..v,(R) ■ -^[1 + sgn(r,, - r)] , (73) 



which goes to zero at large r (see Note roll ). The correct (taking care of the missing factor 
of two in Ref . 



151 ) /3 correction is 
A/f_,,(r,R) = -[i?^(R) - Eo]^ E 5<x,.,5.',., I 



A A— La ^ J / 



rffir 1 




4vr |r-rij 




^ij + ''" |''"ij 


-r| 



Note that also {AI^^,(r, R)) goes to zero at large r. This particular HFM measure needs to 
be shifted g^,ji{r) = g^™{r) + 1. We chose to do the shift as follows: g„,ji{r) = g^™{r) + 
5'™^^(L/2) —g^™{L/2). Nonetheless it is expected to give better results for the on-top value 
of the RDF where the histogram estimator of Eq. (jUj), after the necessary discretization 
of the Dirac delta function, leads, in the measure, to a statistical average divided by zero. 
Moreover it does not suffer from any discretization error and can be calculated for any value 
of r. 

In Fig. |5] we show a comparison for the RDF of the A^ = 162, (^ = 0,/i = oo,rs = 10 
system, calculated in DMC SJ with various kinds of measures. The length of the run was 
always the same 50 blocks of 500 steps each. From the figure one can see that with our 
choice of the /3 correction the HFM measure has the correct average value (coinciding with 
the usual histogram estimator). From the figure it is also evident that the HFM measure is 
much less efficient than the other measures (clearly with a sufficient number of blocks the 
statistical error on the HFM measure can be made small at will). 
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FIG. 5. Shows the RDF of the N = 162, ^ = 0, /x = oo, r^ = 10 system, calculated in DMC SJ with 
various kinds of measures: mixed histogram (mixed), extrapolated histogram (extrapolated), and 
HFM (HFM) with the choice of Eq. (172]) . 

This inefficiency is entirely due to the ZB correction (essential in the DMC calculation). 
From its definition (see Eq. ( 1741) ) one can see that it is the small difference of two large 
terms involving the (extensive) total energy . So the statistical error on the HFM measure is 
completely dominated by that of the /3 part, the a part having statistical errors comparable 
with the ones of the usual histogram estimator, as shown in the left panel of Fig. [61 

E. Choice of the Jastrow 



We noticed that at small r^,/!, and r the variational measure for the unlike RDF, with 
the chosen Jastrow J'l of Eq. ( 136|) . deviates strongly from the mixed one. This is no longer 
so with the modified Jastrow J'2 of Eq. f p9|) which at small /i gives also better variational 
energies (but not for /i — )■ 00 where J7i is better. Note that the Jastrow factor does not 
change the nodes of the wave-function so the energies calculated from the diffusion with 
J'l or J'2 coincide). The extrapolated measures do not change appreciably in the two cases 
apart from near r = 0. In Fig. [6] we show the difference for the two calculations with J7i 
and J2 for the C = 0, r^ = 1, /x = 1 model. From the inset in the left panel we can see that 
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among the two extrapolated measures there is a difference of the order of 0.005. 
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FIG. 6. Unlike RDF for the unpolarized fluid of Pazianii^ at r^ = 1 and f-i = 1 with N = 162. 
On the left the calculation with the Jastrow J'l of Eq. ()36p with various measures: variational 
histogram (variational) and variational HFMv (HFMv) using the estimator of Eq. (j73|) . mixed 
histogram (mixed) and HFM (HFM), and extrapolated histogram (extrapolated). On the right 
the calculation with the Jastrow J'2 of Eq. (j39p with the histogram variational (variational) , mixed 
(mixed), and extrapolated (extrapolated) measures. In the inset is shown the difference between 
the histogram extrapolated measure of the calculation with J'l and the histogram extrapolated 
measure of the calculation with J'2- 10^ Monte Carlo steps were used in the simulations. 



Our results with the two Jastrow factors show that J7i is better than 1/2 for the near- 
Jellium systems (/i large) while J^2 is better than J'l for the near-ideal systems (/i small). 

F. The histogram estimator 

In Fig. [7] we show the DMC results for the histogram extrapolated measure of the RDF 
of our fluid model at C = 0. The time step, r, and number of walkers, n^, were chosen 
according to the indications given in subsection IVIII Al Fig. |8]is for the C = 1 case. 
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FIG. 7. Shows the histogram extrapolated measure for the RDF of a system of 162 unpolarized 
(C = 0) particles calculated using the SJ trial wave-function. The VMC calculation was made of 
10^ steps while the DMC by 10^. The trial wave-function used was of the SJ type with the Jastrow 
Jx of Eq. 
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FIG. 8. Shows the histogram extrapolated measure for the RDF of a system of 147 fully polarized 
(C = 1) particles calculated using the SJ trial wave-function. The VMC calculation was made of 
10^ steps while the DMC by 10^. The trial wave-function used was of the SJ type with the Jastrow 
Ji of Eq. 



In Table |IV] we show the on-top values for the unlike RDF, (7+„(0), of the unpolarized 
system, calculated with the histogram variational, the histogram mixed, the histogram ex- 
trapolated measure, the HFM measure, and the HFM extrapolated measure (of Eq. flM|) ). 
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TABLE IV. Contact values for the unlike RDF of the unpolarized fluid of Paziani^^, at various 
Vg and ^t, from the histogram variational (variational), mixed (mixed), and extrapolated (extrap- 
olated) measures, and the HFM (HFM) and HFM extrapolated (HFM-ext) measures. The trial 
wave- function used was of the SJ type with the Jastrow Ji of Eq. (j36p . The last column gives the 



error cjav = y/cr'^JC/M {(j'^ is the variance, /C the correlation time of the random walk, and M the 
number of Monte Carlo steps) on the HFM measure. 162 particles were used with 10^ x n^ Monte 
Carlo steps. 



Ta 
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variational 


mixed 


extrapolated 
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HFM-ext 


(Jav on HFM 


10 


1/2 


1.085(8) 


1.000(4) 


0.91(1) 


1.0006 


0.9222 


0.03 


10 


1 


0.706(6) 


0.644(3) 


0.582(8) 


0.6474 


0.5949 


0.03 


10 


2 


0.219(4) 


0.182(1) 


0.146(4) 


0.1798 


0.1450 


0.06 


10 


4 


0.053(2) 


0.0506(8) 


0.048(2) 


0.0460 


0.0394 


0.07 


10 


oo 


0.0074(6) 


0.0096(3) 


0.0118(8) 


0.0045 


0.0029 


0.09 


5 


1/2 


1.129(8) 


1.034(3) 


0.94(1) 


1.0277 


0.9381 


0.03 


5 


1 


0.850(7) 


0.796(3) 


0.743(9) 


0.7912 


0.7325 


0.02 


5 


2 


0.448(5) 


0.405(2) 


0.362(6) 


0.4022 


0.3565 


0.02 


5 


4 


0.214(3) 


0.199(1) 


0.184(4) 
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0.1782 


0.03 


5 


oo 


0.080(2) 
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0.080(2) 
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0.0557 
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2 


1/2 


1.158(8) 


1.0618(4) 
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1.0545 
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0.04 


2 


1 


1.003(8) 


0.927(3) 


0.852(9) 
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0.8561 


0.03 
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2 


0.754(7) 


0.697(3) 


0.639(9) 


0.6919 


0.6299 


0.02 


2 
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0.549(6) 


0.511(2) 


0.473(7) 


0.5127 


0.4687 


0.02 


2 


oo 


0.376(4) 


0.349(2) 


0.323(5) 


0.3236 


0.3030 


0.02 
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1.171(8) 


1.077(3) 


0.98(1) 


1.0705 


0.9683 


0.02 




1 


1.077(8) 


0.994(3) 


0.91(1) 


0.9938 


0.9070 


0.02 




2 


0.924(8) 


0.855(3) 


0.787(9) 


0.8640 


0.8053 


0.02 




4 


0.784(7) 


0.730(2) 


0.676(8) 


0.7295 


0.6628 


0.01 




oo 


0.645(6) 


0.602(2) 


0.560(7) 


0.5771 


0.5263 


0.01 



IX. CONCLUSIONS 



We studied through Variational and Diffusion Monte Carlo techniques the fluid of spin 
one-half particles interacting with the bare pair-potential Vfj,{r) = eii{^r)/r and immersed 
in a uniform counteracting background. When // — ;■ oo the system reduces to the Jellium 
model whereas when // — )■ it reduces to the ideal Fermi gas. We performed a detailed 
analysis of the spin-resolved Radial Distribution Function for this system as a function of 
the density parameter r^ = 1, 2, 5, 10 and the penetrability parameter // = 1/2, 1, 2, 4, oo at 
two values of the polarization, ^ = 0, 1. 
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Initially we carefully fine tuned our DMC calculation determining the optimal values for 
the time step r and the number of walkers n^ for each value of the density parameter r^. 
Increasing the system size A^ the RDF extends its range [0, Tmax] at larger r^ax- We estimated 
that for A^ > 66 the size dependence of the RDF is lower than 2%. As a compromise 
between computational cost and reduction of the size effects, the largest uncontrolled source 
of uncertainty on our RDF measurements, we chose to perform the RDF calculation with 
A^ = 162 in the unpolarized case and A^ = 147 in the polarized case. 

We calculated the RDF using two different routes: through the usual histogram estimator 
and through a particular HFM measure. As expected, in the VMC calculations the HFMv 
estimator gives better results for the on-top value of the RDF. In the DMC calculation the 
inclusion of the /3 correction (which must be omitted in the VMC calculation) is indispens- 
able. Moreover the ZV estimator is zero for r > Vmax so it has to be shifted by +1. From 
our variational and fixed nodes diffusion Monte Carlo experiments turns out that although 
in the variational measure the average of the histogram estimator agrees with the average of 



the HFMv estimator within the square root of the variance of the average o"av = \/o^C/Af, 
where a^ is the variance, /C the correlation time of the random walk, and A/" the number of 
Monte Carlo steps, and the two (Xav are comparable, in the diffusion experiment, where one 
has to add the /3 correction not to bias the average, the Hellmann and Feynman measure 
has an average in agreement with the one of the histogram estimator but the a^v increases. 
This is to be expected from the extensive nature of the /3 correction in which the energy 
appears. Of course the averages from the extrapolated Hellmann and Feynman measure and 
the extrapolated measure for the histogram estimator also agree. 

In the simulation, for the Coulomb case, fi — > oo, we made extrapolations in time step 
and number of walkers for each value of r^. Given a relative precision S^o = Aeo/Cp, where 
Co = {El) f/N, Aeo is the statistical error on cq, and e^ is the exchange energy, we set 
as our target relative precision 5^^ = 10~^%. The extrapolated values of the time step 
and number of walkers was then used for all other values of fi. We chose the trial wave 
function of the Bijl-Dingle-Jastrow^ form as a product of Slater determinats and a Jastrow 
factor. The pseudo potential was chosen as in Ref.— , 1/2, which is expected to give better 
results for Jellium. Comparison with the simulation of the unpolarized fluid at r^ = 1 and 
/i = 1 with the pseudo potential of Ref.—, jTi, for which the trial wave function becomes 
the exact ground state wave function in the /x — )■ limit, show that the two extrapolated 
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measures of the unlike histogram estimator differ one from the other by less than 7 x 10~ , 
the largest difference being at contact (see Fig. 1). The use of more sophisticated trial wave 
functions, taking into account the effect of backflow and three-body correlations, is found to 
affect the measure by even less in the range of densities considered. For the same reason we 
discarded the use of the twist-averaged boundary conditions^ and only worked with periodic 
boundary conditions. In Table IIVI we compare the contact values of the unlike RDF of the 
unpolarized fluid at various r^ and fi from the measures of the histogram estimator and the 
HFM measures. We see that there is disagreement between the measure from the histogram 
estimator and the HFM measure only in the Coulomb /i — )■ oo case at r^ = 1,2. 

Our results complement the ones of Paziani et ai.— which only reported a limited number 
of RDF data. We plan, in the future, to complete the calculation at intermediate polar- 
izations, < C < 1, complementing the work of Ortiz and Ballone^, and Ceperley and 
coworkers^^. 

We believe it is still an open problem the one of determining the relationship between 
the choice of the auxiliary function, the bias it introduces in the Hellmann and Feynman 
measure, and the variance of this measure. 

Appendix A: Jastrow, backflow, and three-body 



In terms of the stochastic process governed by /(R, t) one can write, using Kac theorent^i; 



63M 



/"rfR/(R,r) = /exp - f dtEL(R') 



(Al) 

DRW 



where (. . .)drw means averaging with respect to the diffusing and drifting random walk. 
Choosing a complete set of orthonormal wave-functions \l/j we can write for the true time 
dependent many-body wave-function 

0(R,r) = ^*i(R) /"rfR'*,(R')0(R',r) ^ ^(R) /"rfR/(R,r) 
= *(R)/exp - f dtELiU! 



(A2) 

DRW 



'0 

where \1/ is the wave-function, of the set, of maximum overlap with the true ground-state, 
the trial wave-function. Assuming that at time zero we are already close to the stationary 
solution, for sufficiently small r we can approximate 

e-^^^(^^) . (A3) 



( exp 



dtELiR') 
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f-*^ 



DRW 



By antisymmetrising we get the Fermion wave-function 

0^(R,r)^^[e-^^-Wv[/(R)] , (A4) 

where given a function /(R) we define the operator (a symmetry of the Hamiltonian) 

^[/(R)] = -^^(-l)^/(PR), (A5) 

^ p 

here Np = N^lN^l is the total number of allowed permutations P. 

This is called the local energy method to improve a trial wave-function. Suppose we start 

from a simple unsymmetrical product of single particle plane waves of A^+ spin-up particles 

with k < kp occupied and A^_ spin-up particles with k < kp occupied, for the zeroth order 

trial wave-function. Equation (IA4p will give us a first order wave-function of the Slater- 

Jastrow type (see equation (!27|) ). If we start from an unsymmetrical Hartree-Jastrow trial 

wave-function the local energy with the Jastrow factor has the form 

2' 



El = V-\Y, 



-kl - 2zki ■ V, ^ M(r,fc) - V,' ^ M(r,fc) + 

j<k j<k 



j<k 



(A6) 



where V = V(R) is the total potential energy and r^ = \rij\ = |rj — rj\. Then the antisym- 
metrized second order wave- function has the form in Eq. (1401) . which includes backflow (see 
the third term), which is the correction inside the determinant and which affects the nodes, 
and three-body boson-like correlations (see last term) which do not affect the nodes. 
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